Analysis of Freesurfer parcellation Genetic Template, surface area and cortical thickness differences

library(easypackages)
libraries("here","heplots","ggplot2","ggseg","ggsegExtra","ggsegChen","tidyverse","patchwork","pheatmap","psych")
source(here("code","cohens_d.R"))
source(here("code","get_ggColorHue.R"))

datapath = here("data","tidy")
plotdir = here("plots")

pheno_data_gex =read.csv(file.path(datapath,"labelData_allGEXMRIsubs.csv"))
pheno_data_mri =read.csv(file.path(datapath,"labelData_all_MRI.csv"))


# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)

fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)

region_names = colnames(fs_data)[5:ncol(fs_data)]
region_names = c("CortexVol","SAtotal","CTmean",region_names)

fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")

new_fs_data = fs_data

colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
                 "tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
                 "pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use

for (icol in 1:length(region_names)){
  
  y_var = region_names[icol]
  
  if(is.element(y_var,c("CortexVol","SAtotal","CTmean"))){
    # make formula
    form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
    full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
      covname2use = c("sexMale","scan_age")
  }else{
    # make formula
    form2use = as.formula(sprintf("%s ~ Group + CortexVol + sex + scan_age", y_var))
    full_model = model.matrix(~0+ Group + CortexVol + sex + scan_age, data = fs_data)
    covname2use = c("CortexVol","sexMale","scan_age")
  }
  
  # run linear model
  mod2use = lm(formula = form2use, data = fs_data)
  
  # run an ANOVA
  res2use = anova(mod2use)
  
  # grab ANOVA results and store into output_res
  output_res[y_var, "Fstat"] = res2use["Group", "F value"]
  output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
  
  # compute partial eta-squared as the effect size of the interaction effect
  eta_sq_res = etasq(mod2use)
  output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
  
  
    # remove covariates before effect size computation
    beta1 = mod2use$coefficients[covname2use, drop = FALSE]
    beta1[is.na(beta1)] = 0

    new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
    
    # compute effect sizes
    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"] 
    
    mask1 = fs_data$Group=="Poor"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 


    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="Poor"
    output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])

    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 

}

output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
output_res[order(output_res$pval),]
##                                         Fstat         pval        etasq
## CortexVol                          18.7684243 3.563168e-08 0.0609961017
## anteromedialtemporal.lh.area       14.0080609 2.095489e-06 0.1018626300
## SAtotal                            12.7918734 6.069546e-06 0.0384031239
## superiorparietal.rh.area           10.6221181 4.214024e-05 0.0406886255
## posterolateraltemporal.rh.area     10.1768804 6.297412e-05 0.0532927353
## inferiorparietal.lh.thickness       9.0672894 1.726457e-04 0.0855251058
## superiorparietal.lh.area            8.8900453 2.030244e-04 0.0469626440
## middletemporal.rh.thickness         7.2377106 9.322992e-04 0.0567858094
## dorsomedialfrontal.lh.area          6.9207431 1.252445e-03 0.0500938981
## perisylvian.lh.thickness            5.6069622 4.299179e-03 0.0280653457
## middletemporal.lh.thickness         5.1161171 6.843551e-03 0.0160809808
## inferiorparietal.rh.thickness       5.1077155 6.898359e-03 0.0452499419
## orbitofrontal.rh.area               4.7199377 9.975975e-03 0.0884379862
## dorsomedialfrontal.rh.area          4.6591244 1.057149e-02 0.0500313212
## medialprefrontal.rh.thickness       4.4135056 1.336577e-02 0.0035100067
## medialprefrontal.lh.thickness       4.3594368 1.407504e-02 0.0049069724
## inferiorparietal.rh.area            4.1241893 1.763166e-02 0.0388947624
## orbitofrontal.lh.area               3.9924278 2.000762e-02 0.0817033669
## dorsolateralprefrotal.lh.thickness  3.7760453 2.463286e-02 0.0010251563
## posterolateraltemporal.lh.area      3.6074681 2.897450e-02 0.0112156872
## occipital.lh.area                   3.5466894 3.072299e-02 0.0234903844
## temporalpole.lh.thickness           3.1731609 4.407572e-02 0.0117222407
## ventromedialoccipital.lh.thickness  2.8410823 6.081898e-02 0.0428891443
## superiorparietal.rh.thickness       2.8250048 6.177622e-02 0.0223598673
## anteromedialtemporal.rh.area        2.6898314 7.045102e-02 0.0075296610
## medialtemporal.lh.thickness         2.5864636 7.790740e-02 0.0314594717
## temporalpole.rh.thickness           2.5859020 7.795001e-02 0.0019890549
## inferiorparietal.lh.area            2.5563368 8.022714e-02 0.0301583009
## CTmean                              2.5204136 8.307086e-02 0.0011589273
## dorsolateralprefrotal.rh.thickness  2.4716397 8.712967e-02 0.0099715010
## ventralfrontal.rh.thickness         2.3301281 1.000286e-01 0.0204949358
## motor_premotor_SMA.rh.thickness     2.1554571 1.186455e-01 0.0549536354
## perisylvian.rh.thickness            2.1521394 1.190311e-01 0.0131780702
## superiortemporal.rh.area            2.0552479 1.308692e-01 0.0056249313
## motor_premotor_SMA.lh.thickness     1.6863034 1.879341e-01 0.0782779395
## dorsolateralprefrontal.rh.area      1.4145993 2.455442e-01 0.0080141342
## occipital.lh.thickness              1.2841529 2.792522e-01 0.0109164163
## occipital.rh.thickness              1.2733925 2.822334e-01 0.0057228488
## precuneus.rh.area                   1.2492498 2.890396e-01 0.0107147043
## motorpremotor.rh.area               1.0954063 3.364843e-01 0.0149046475
## ventromedialoccipital.rh.thickness  1.0815033 3.411419e-01 0.0157302561
## parsopercularis.lh.area             1.0271582 3.599817e-01 0.0107116670
## ventralfrontal.lh.thickness         0.9478539 3.893778e-01 0.0136577437
## dorsolateralprefrontal.lh.area      0.5413956 5.828213e-01 0.0153294284
## precuneus.lh.area                   0.4888912 6.140674e-01 0.0038323834
## medialtemporal.rh.thickness         0.2924523 7.467628e-01 0.0131748655
## motorpremotor.lh.area               0.2654045 7.671765e-01 0.0029461200
## occipital.rh.area                   0.2082300 8.122034e-01 0.0001698168
## parsopercularis.rh.area             0.1699284 8.438520e-01 0.0001588197
## superiortemporal.lh.area            0.0788472 9.242110e-01 0.0020092765
## superiorparietal.lh.thickness       0.0213230 9.789050e-01 0.0004093216
##                                            fdrq es.TD_vs_Good es.TD_vs_Poor
## CortexVol                          1.817216e-06  -0.168185504  -0.627477606
## anteromedialtemporal.lh.area       5.343498e-05   0.073908193  -0.716130093
## SAtotal                            1.031823e-04  -0.086079474  -0.471191290
## superiorparietal.rh.area           5.372881e-04  -0.239268265   0.266474769
## posterolateraltemporal.rh.area     6.423360e-04   0.166221696  -0.413272621
## inferiorparietal.lh.thickness      1.467488e-03   0.402056404  -0.376038602
## superiorparietal.lh.area           1.479178e-03  -0.338833602   0.223427618
## middletemporal.rh.thickness        5.943408e-03  -0.602325506  -0.335190443
## dorsomedialfrontal.lh.area         7.097189e-03  -0.413580399   0.117701311
## perisylvian.lh.thickness           2.192581e-02  -0.394967263  -0.059355222
## middletemporal.lh.thickness        2.931803e-02  -0.078034977  -0.363741804
## inferiorparietal.rh.thickness      2.931803e-02   0.438813968  -0.073227995
## orbitofrontal.rh.area              3.851042e-02   0.703440285   0.689021478
## dorsomedialfrontal.rh.area         3.851042e-02  -0.011856131   0.560176404
## medialprefrontal.rh.thickness      4.486419e-02   0.079704943   0.150731011
## medialprefrontal.lh.thickness      4.486419e-02  -0.005454202   0.155007839
## inferiorparietal.rh.area           5.289497e-02  -0.496123250  -0.002619246
## orbitofrontal.lh.area              5.668827e-02   0.313080777   0.746115082
## dorsolateralprefrotal.lh.thickness 6.611979e-02   0.056890925   0.081060057
## posterolateraltemporal.lh.area     7.388497e-02   0.074128125  -0.187712338
## occipital.lh.area                  7.461298e-02   0.395346960   0.163865789
## temporalpole.lh.thickness          1.021755e-01  -0.011829570  -0.257701161
## ventromedialoccipital.lh.thickness 1.312745e-01  -0.415214100  -0.555129522
## superiorparietal.rh.thickness      1.312745e-01   0.308387021  -0.015645109
## anteromedialtemporal.rh.area       1.437201e-01  -0.067557401  -0.230844671
## medialtemporal.lh.thickness        1.460901e-01   0.100526613   0.450751114
## temporalpole.rh.thickness          1.460901e-01  -0.102905181  -0.020598602
## inferiorparietal.lh.area           1.460901e-01  -0.297574970  -0.448702546
## CTmean                             1.460901e-01   0.022296408   0.086644973
## dorsolateralprefrotal.rh.thickness 1.481204e-01   0.140067779  -0.103150861
## ventralfrontal.rh.thickness        1.645632e-01  -0.323190343  -0.295450988
## motor_premotor_SMA.rh.thickness    1.839571e-01   0.313202751   0.608448832
## perisylvian.rh.thickness           1.839571e-01  -0.119473742  -0.292404471
## superiortemporal.rh.area           1.963038e-01   0.182570796   0.138964153
## motor_premotor_SMA.lh.thickness    2.738468e-01   0.573830868   0.746481899
## dorsolateralprefrontal.rh.area     3.478543e-01   0.167982081  -0.040570128
## occipital.lh.thickness             3.779749e-01  -0.111820856   0.148339120
## occipital.rh.thickness             3.779749e-01  -0.013600907   0.166091886
## precuneus.rh.area                  3.779749e-01  -0.117983093   0.134026057
## motorpremotor.rh.area              4.243472e-01  -0.207271672  -0.303815461
## ventromedialoccipital.rh.thickness 4.243472e-01  -0.302852519  -0.232545925
## parsopercularis.lh.area            4.371207e-01   0.026319646  -0.249586445
## ventralfrontal.lh.thickness        4.618202e-01  -0.264409577  -0.259488718
## dorsolateralprefrontal.lh.area     6.755429e-01   0.171882416   0.316840463
## precuneus.lh.area                  6.959431e-01  -0.077260554   0.076319613
## medialtemporal.rh.thickness        8.279326e-01   0.089540153   0.324278245
## motorpremotor.lh.area              8.324682e-01  -0.010298772   0.112012146
## occipital.rh.area                  8.629661e-01   0.030998595   0.028420544
## parsopercularis.rh.area            8.782949e-01  -0.027358771  -0.026814437
## superiortemporal.lh.area           9.426953e-01  -0.114992084  -0.083632069
## superiorparietal.lh.thickness      9.789050e-01   0.023147126   0.057755463
##                                    es.Good_vs_Poor tstat.TD_vs_Good
## CortexVol                             -0.416871265      -0.92763624
## anteromedialtemporal.lh.area          -0.758745394       0.37243309
## SAtotal                               -0.372984599      -0.53069497
## superiorparietal.rh.area               0.547055816      -1.56200620
## posterolateraltemporal.rh.area        -0.604517309       1.15445234
## inferiorparietal.lh.thickness         -0.744128402       2.52598615
## superiorparietal.lh.area               0.528719718      -2.00216978
## middletemporal.rh.thickness            0.292282835      -3.22118075
## dorsomedialfrontal.lh.area             0.538454727      -2.36785475
## perisylvian.lh.thickness               0.344432195      -2.31179712
## middletemporal.lh.thickness           -0.226876400      -0.59981901
## inferiorparietal.rh.thickness         -0.523112641       2.46021689
## orbitofrontal.rh.area                 -0.022666089       3.81615464
## dorsomedialfrontal.rh.area             0.481823664      -0.03266780
## medialprefrontal.rh.thickness          0.075467838       0.62367574
## medialprefrontal.lh.thickness          0.148087974      -0.10676434
## inferiorparietal.rh.area               0.396157491      -2.97475639
## orbitofrontal.lh.area                  0.472178056       1.54212952
## dorsolateralprefrotal.lh.thickness     0.021731886       0.08491989
## posterolateraltemporal.lh.area        -0.267730933       0.36341567
## occipital.lh.area                     -0.219531034       2.06315786
## temporalpole.lh.thickness             -0.229334676      -0.05302691
## ventromedialoccipital.lh.thickness    -0.054958989      -2.36766848
## superiorparietal.rh.thickness         -0.350452936       1.98361558
## anteromedialtemporal.rh.area          -0.144985360      -0.48454790
## medialtemporal.lh.thickness            0.348994398       0.52586747
## temporalpole.rh.thickness              0.082064561      -0.68884090
## inferiorparietal.lh.area              -0.177660163      -1.59491429
## CTmean                                 0.059512817       0.06023806
## dorsolateralprefrotal.rh.thickness    -0.269073157       0.28318744
## ventralfrontal.rh.thickness            0.049662712      -2.02068983
## motor_premotor_SMA.rh.thickness        0.317960692       1.78735887
## perisylvian.rh.thickness              -0.180757516      -0.92800055
## superiortemporal.rh.area              -0.041254182       1.26918573
## motor_premotor_SMA.lh.thickness        0.132796743       3.26084740
## dorsolateralprefrontal.rh.area        -0.225924020       1.12326898
## occipital.lh.thickness                 0.273712261      -0.76184410
## occipital.rh.thickness                 0.172930032       0.02467738
## precuneus.rh.area                      0.288530088      -0.69967304
## motorpremotor.rh.area                 -0.101723418      -1.08822596
## ventromedialoccipital.rh.thickness     0.095756438      -1.51414221
## parsopercularis.lh.area               -0.235320467       0.06273803
## ventralfrontal.lh.thickness           -0.008185191      -1.58638447
## dorsolateralprefrontal.lh.area         0.153830025       1.12998643
## precuneus.lh.area                      0.161333323      -0.29115018
## medialtemporal.rh.thickness            0.193392174       0.62787481
## motorpremotor.lh.area                  0.122913072      -0.03900672
## occipital.rh.area                      0.001588241      -0.08532334
## parsopercularis.rh.area                0.003390924       0.01719381
## superiortemporal.lh.area               0.030124858      -0.30438030
## superiorparietal.lh.thickness          0.026537911       0.53721444
##                                    tstat.TD_vs_Poor tstat.Good_vs_Poor
## CortexVol                              -3.059262385         2.51015442
## anteromedialtemporal.lh.area           -3.515199051         4.06483461
## SAtotal                                -2.131564929         2.19003479
## superiorparietal.rh.area                1.984471136        -2.77322638
## posterolateraltemporal.rh.area         -2.047477385         3.54666956
## inferiorparietal.lh.thickness          -2.131633301         4.08238560
## superiorparietal.lh.area                1.551859981        -2.78194967
## middletemporal.rh.thickness            -1.778064339        -1.60535282
## dorsomedialfrontal.lh.area              0.509988831        -3.41232980
## perisylvian.lh.thickness               -0.296281168        -2.01527339
## middletemporal.lh.thickness            -1.040308667         1.32359190
## inferiorparietal.rh.thickness          -0.867470744         2.38090177
## orbitofrontal.rh.area                   3.758467841         0.06798017
## dorsomedialfrontal.rh.area              3.017287570        -2.63515835
## medialprefrontal.rh.thickness           0.871125758        -0.31661320
## medialprefrontal.lh.thickness           0.819719437        -0.78027216
## inferiorparietal.rh.area                0.499151365        -1.80659810
## orbitofrontal.lh.area                   3.924457099        -3.01607583
## dorsolateralprefrotal.lh.thickness      0.722697259        -0.12019995
## posterolateraltemporal.lh.area         -0.612513896         1.73587140
## occipital.lh.area                       0.228185483         0.88449841
## temporalpole.lh.thickness              -1.353811398         1.51173501
## ventromedialoccipital.lh.thickness     -2.908615860        -0.02696405
## superiorparietal.rh.thickness          -0.149573963         1.97211361
## anteromedialtemporal.rh.area           -0.893081728         0.90645129
## medialtemporal.lh.thickness             2.275557563        -1.70721752
## temporalpole.rh.thickness              -0.065063096        -0.25861361
## inferiorparietal.lh.area               -2.228747825         1.12256572
## CTmean                                  0.003809356        -0.36930388
## dorsolateralprefrotal.rh.thickness     -0.083335625         1.44315454
## ventralfrontal.rh.thickness            -1.470130934        -0.01116362
## motor_premotor_SMA.rh.thickness         3.062774577        -1.85472165
## perisylvian.rh.thickness               -1.210572588         0.71872792
## superiortemporal.rh.area                0.356212566         0.34908797
## motor_premotor_SMA.lh.thickness         3.933261727        -0.91827562
## dorsolateralprefrontal.rh.area         -0.331176983         1.28323113
## occipital.lh.thickness                  0.315153460        -1.87591783
## occipital.rh.thickness                  0.828635994        -0.78044635
## precuneus.rh.area                       1.160264079        -1.74582006
## motorpremotor.rh.area                  -1.779218478         0.37142742
## ventromedialoccipital.rh.thickness     -1.134943163        -0.54892205
## parsopercularis.lh.area                -1.564937351         1.22474572
## ventralfrontal.lh.thickness            -1.209652541         0.50358676
## dorsolateralprefrontal.lh.area          1.193483568        -0.85504869
## precuneus.lh.area                       0.844134840        -0.70718029
## medialtemporal.rh.thickness             1.638318531        -0.72499110
## motorpremotor.lh.area                   0.837840504        -0.34973418
## occipital.rh.area                      -0.333265662        -0.45488366
## parsopercularis.rh.area                -0.860511040        -0.10119179
## superiortemporal.lh.area               -0.762767692        -0.15225809
## superiorparietal.lh.thickness          -0.143367881        -0.20504914
##                                    pval.TD_vs_Good pval.TD_vs_Poor
## CortexVol                             0.3553559051    0.0026877095
## anteromedialtemporal.lh.area          0.7101961158    0.0006040535
## SAtotal                               0.5965569661    0.0348955186
## superiorparietal.rh.area              0.1207947235    0.0492929914
## posterolateraltemporal.rh.area        0.2504996384    0.0426100630
## inferiorparietal.lh.thickness         0.0127767912    0.0349040119
## superiorparietal.lh.area              0.0474136371    0.1231083165
## middletemporal.rh.thickness           0.0016246125    0.0777132171
## dorsomedialfrontal.lh.area            0.0194129637    0.6109174538
## perisylvian.lh.thickness              0.0224120292    0.7674845535
## middletemporal.lh.thickness           0.5497040259    0.3001127370
## inferiorparietal.rh.thickness         0.0152403869    0.3872703845
## orbitofrontal.rh.area                 0.0002113866    0.0002564376
## dorsomedialfrontal.rh.area            0.9739911935    0.0030651648
## medialprefrontal.rh.thickness         0.5339679967    0.3852791198
## medialprefrontal.lh.thickness         0.9151457399    0.4138649853
## inferiorparietal.rh.area              0.0035153952    0.6185099704
## orbitofrontal.lh.area                 0.1255499164    0.0001396919
## dorsolateralprefrotal.lh.thickness    0.9324598646    0.4711541563
## posterolateraltemporal.lh.area        0.7169034100    0.5412595960
## occipital.lh.area                     0.0411514717    0.8198578416
## temporalpole.lh.thickness             0.9577944449    0.1781276166
## ventromedialoccipital.lh.thickness    0.0194223111    0.0042658904
## superiorparietal.rh.thickness         0.0494727774    0.8813307850
## anteromedialtemporal.rh.area          0.6288383274    0.3734510978
## medialtemporal.lh.thickness           0.5999046095    0.0244961519
## temporalpole.rh.thickness             0.4921898421    0.9482229469
## inferiorparietal.lh.area              0.1132367114    0.0275351297
## CTmean                                0.9520607209    0.9969663319
## dorsolateralprefrotal.rh.thickness    0.7774978971    0.9337118827
## ventralfrontal.rh.thickness           0.0454314672    0.1439240155
## motor_premotor_SMA.rh.thickness       0.0762835306    0.0026621984
## perisylvian.rh.thickness              0.3551815765    0.2282388836
## superiortemporal.rh.area              0.2067142418    0.7222544192
## motor_premotor_SMA.lh.thickness       0.0014289517    0.0001351941
## dorsolateralprefrontal.rh.area        0.2634589483    0.7410399734
## occipital.lh.thickness                0.4475767309    0.7531466001
## occipital.rh.thickness                0.9803513384    0.4088176208
## precuneus.rh.area                     0.4854212466    0.2480514687
## motorpremotor.rh.area                 0.2785725200    0.0775229375
## ventromedialoccipital.rh.thickness    0.1324942477    0.2584714928
## parsopercularis.lh.area               0.9500745288    0.1200101061
## ventralfrontal.lh.thickness           0.1151586014    0.2285906765
## dorsolateralprefrontal.lh.area        0.2606284515    0.2348369502
## precuneus.lh.area                     0.7714151924    0.4001329423
## medialtemporal.rh.thickness           0.5312222760    0.1037549440
## motorpremotor.lh.area                 0.9689467714    0.4036463233
## occipital.rh.area                     0.9321397658    0.7394665581
## parsopercularis.rh.area               0.9863091943    0.3910795681
## superiortemporal.lh.area              0.7613403104    0.4469732857
## superiorparietal.lh.thickness         0.5920669333    0.8862199452
##                                    pval.Good_vs_Poor fdrq.TD_vs_Good
## CortexVol                               1.334511e-02      0.69704428
## anteromedialtemporal.lh.area            8.475100e-05      0.96215984
## SAtotal                                 3.037760e-02      0.87414672
## superiorparietal.rh.area                6.408005e-03      0.35564245
## posterolateraltemporal.rh.area          5.514425e-04      0.58419158
## inferiorparietal.lh.thickness           7.930211e-05      0.12381723
## superiorparietal.lh.area                6.248267e-03      0.19408551
## middletemporal.rh.thickness             1.109600e-01      0.02761841
## dorsomedialfrontal.lh.area              8.703711e-04      0.12381723
## perisylvian.lh.thickness                4.603818e-02      0.12700150
## middletemporal.lh.thickness             1.880737e-01      0.87414672
## inferiorparietal.rh.thickness           1.879256e-02      0.12381723
## orbitofrontal.rh.area                   9.459109e-01      0.01078072
## dorsomedialfrontal.rh.area              9.482980e-03      0.98630919
## medialprefrontal.rh.thickness           7.520695e-01      0.87414672
## medialprefrontal.lh.thickness           4.367176e-01      0.98630919
## inferiorparietal.rh.area                7.325004e-02      0.04482129
## orbitofrontal.lh.area                   3.107251e-03      0.35564245
## dorsolateralprefrotal.lh.thickness      9.045193e-01      0.98630919
## posterolateraltemporal.lh.area          8.507058e-02      0.96215984
## occipital.lh.area                       3.781391e-01      0.19408551
## temporalpole.lh.thickness               1.331457e-01      0.98630919
## ventromedialoccipital.lh.thickness      9.785318e-01      0.12381723
## superiorparietal.rh.thickness           5.082244e-02      0.19408551
## anteromedialtemporal.rh.area            3.664551e-01      0.89085430
## medialtemporal.lh.thickness             9.028439e-02      0.87414672
## temporalpole.rh.thickness               7.963622e-01      0.86557524
## inferiorparietal.lh.area                2.637909e-01      0.35564245
## CTmean                                  7.125261e-01      0.98630919
## dorsolateralprefrotal.rh.thickness      1.514984e-01      0.96713153
## ventralfrontal.rh.thickness             9.911108e-01      0.19408551
## motor_premotor_SMA.rh.thickness         6.601136e-02      0.27789000
## perisylvian.rh.thickness                4.736608e-01      0.69704428
## superiortemporal.rh.area                7.276154e-01      0.52712132
## motor_premotor_SMA.lh.thickness         3.602573e-01      0.02761841
## dorsolateralprefrontal.rh.area          2.018042e-01      0.58419158
## occipital.lh.thickness                  6.301777e-02      0.84542271
## occipital.rh.thickness                  4.366155e-01      0.98630919
## precuneus.rh.area                       8.331901e-02      0.86557524
## motorpremotor.rh.area                   7.109531e-01      0.59196660
## ventromedialoccipital.rh.thickness      5.840460e-01      0.35564245
## parsopercularis.lh.area                 2.229926e-01      0.98630919
## ventralfrontal.lh.thickness             6.154453e-01      0.35564245
## dorsolateralprefrontal.lh.area          3.941733e-01      0.58419158
## precuneus.lh.area                       4.807814e-01      0.96713153
## medialtemporal.rh.thickness             4.698234e-01      0.87414672
## motorpremotor.lh.area                   7.271315e-01      0.98630919
## occipital.rh.area                       6.499883e-01      0.98630919
## parsopercularis.rh.area                 9.195617e-01      0.98630919
## superiortemporal.lh.area                8.792308e-01      0.96713153
## superiorparietal.lh.thickness           8.378699e-01      0.87414672
##                                    fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## CortexVol                              0.022331915       0.075622317
## anteromedialtemporal.lh.area           0.007701682       0.002161151
## SAtotal                                0.148342050       0.140841621
## superiorparietal.rh.area               0.179567326       0.046686893
## posterolateraltemporal.rh.area         0.167162555       0.009374523
## inferiorparietal.lh.thickness          0.148342050       0.002161151
## superiorparietal.lh.area               0.330448639       0.046686893
## middletemporal.rh.thickness            0.247710879       0.282947979
## dorsomedialfrontal.lh.area             0.788600212       0.011097232
## perisylvian.lh.thickness               0.869815827       0.195662284
## middletemporal.lh.thickness            0.566879614       0.417033003
## inferiorparietal.rh.thickness          0.603060407       0.095842051
## orbitofrontal.rh.area                  0.004359440       0.984519528
## dorsomedialfrontal.rh.area             0.022331915       0.060453999
## medialprefrontal.rh.thickness          0.603060407       0.891989414
## medialprefrontal.lh.thickness          0.603060407       0.696018673
## inferiorparietal.rh.area               0.788600212       0.233484491
## orbitofrontal.lh.area                  0.003562144       0.031693961
## dorsolateralprefrotal.lh.thickness     0.649428702       0.977034258
## posterolateraltemporal.lh.area         0.726427353       0.241033306
## occipital.lh.area                      0.908972824       0.665003166
## temporalpole.lh.thickness              0.432595640       0.323353722
## ventromedialoccipital.lh.thickness     0.027195051       0.991110843
## superiorparietal.rh.thickness          0.941608692       0.199380340
## anteromedialtemporal.rh.area           0.603060407       0.665003166
## medialtemporal.lh.thickness            0.138811527       0.242342319
## temporalpole.rh.thickness              0.967187406       0.923056146
## inferiorparietal.lh.area               0.140429161       0.517435932
## CTmean                                 0.996966332       0.883532933
## dorsolateralprefrotal.rh.thickness     0.967187406       0.351200902
## ventralfrontal.rh.thickness            0.367006240       0.991110843
## motor_premotor_SMA.rh.thickness        0.022331915       0.224438630
## perisylvian.rh.thickness               0.499028519       0.700567122
## superiortemporal.rh.area               0.869815827       0.883532933
## motor_premotor_SMA.lh.thickness        0.003562144       0.665003166
## dorsolateralprefrontal.rh.area         0.869815827       0.428833975
## occipital.lh.thickness                 0.869815827       0.224438630
## occipital.rh.thickness                 0.603060407       0.696018673
## precuneus.rh.area                      0.506024996       0.241033306
## motorpremotor.rh.area                  0.247710879       0.883532933
## ventromedialoccipital.rh.thickness     0.507001774       0.827398532
## parsopercularis.lh.area                0.330448639       0.454904910
## ventralfrontal.lh.thickness            0.499028519       0.848316434
## dorsolateralprefrontal.lh.area         0.499028519       0.670094525
## precuneus.lh.area                      0.603060407       0.700567122
## medialtemporal.rh.thickness            0.311264832       0.700567122
## motorpremotor.lh.area                  0.603060407       0.883532933
## occipital.rh.area                      0.869815827       0.872352757
## parsopercularis.rh.area                0.603060407       0.977034258
## superiortemporal.lh.area               0.633212155       0.974799407
## superiorparietal.lh.thickness          0.941608692       0.949585942
regions2use = colnames(fs_data)[5:(ncol(fs_data)-4)]
full_regnames = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal")

sa_ordering = c(1,12,7,10,
                4,6,9,3,
                8,11,2,5)
ct_ordering = c(1,7,6,3,
                10,4,11,5,
                8,2,12,9)
gensim_order = c(sa_ordering,sa_ordering, ct_ordering,ct_ordering)

df4plot = output_res[regions2use,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")]
df4plot$phenotype = "CT"
df4plot$phenotype[1:24] = "SA"
df4plot$phenotype = factor(df4plot$phenotype) 

df4plot$hemisphere = "RH"
df4plot$hemisphere[1:12] = "LH"
df4plot$hemisphere[25:36] = "LH"
df4plot$hemisphere = factor(df4plot$hemisphere) 

df4plot$full_region_names = factor(full_regnames)
df4plot$gensim_gradient = gensim_order

ct_df = subset(df4plot, df4plot$phenotype=="CT")
sa_df = subset(df4plot, df4plot$phenotype=="SA")

ct_df$cluster = "dorsal"
ct_df$cluster[c(6,7,8,9,10,11,12)] = "ventral"
ct_df$cluster = factor(ct_df$cluster)

sa_df$cluster = "anterior"
sa_df$cluster[c(6,7,8,9,10,11,12)] = "posterior"
sa_df$cluster = factor(sa_df$cluster)

dotSize = 10
fontSize = 20

p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Good
## t = -1.838, df = 22, p-value = 0.07961
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.66968163  0.04520194
## sample estimates:
##        cor 
## -0.3648475
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDGood_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 0.55381, df = 20.769, p-value = 0.5856
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1553724  0.2680582
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           0.007270414          -0.049072461
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Poor
## t = -0.88864, df = 22, p-value = 0.3838
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5483668  0.2348869
## sample estimates:
##        cor 
## -0.1861478
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_TD_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = -1.1479, df = 18.319, p-value = 0.2658
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3867892  0.1132490
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           -0.05450477            0.08226532
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

cor.test(ct_df$gensim_gradient, ct_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.Good_vs_Poor
## t = 0.80097, df = 22, p-value = 0.4317
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2521901  0.5353739
## sample estimates:
##       cor 
## 0.1683315
p = ggplot(data = ct_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"dvGradient_ASDGood_vs_ASDPoor_CT_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = -1.7772, df = 15.521, p-value = 0.09513
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.41512987  0.03702278
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           -0.05560184            0.13345171
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Good
## t = -0.84267, df = 22, p-value = 0.4085
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.541591  0.243973
## sample estimates:
##        cor 
## -0.1768267
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDGood_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Good ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 1.1527, df = 15.531, p-value = 0.2665
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1017404  0.3429604
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.03042461             -0.09018539
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Poor
## t = -0.82792, df = 22, p-value = 0.4166
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5393989  0.2468832
## sample estimates:
##       cor 
## -0.173825
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_TD_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.TD_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = 1.6011, df = 10.357, p-value = 0.1394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09637807  0.59688583
## sample estimates:
##  mean in group anterior mean in group posterior 
##               0.1092268              -0.1410271
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"geneSimGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

cor.test(sa_df$gensim_gradient, sa_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.Good_vs_Poor
## t = -0.13004, df = 22, p-value = 0.8977
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4263442  0.3799303
## sample estimates:
##         cor 
## -0.02771394
p = ggplot(data = sa_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
ggsave(filename = file.path(plotdir,"apGradient_ASDGood_vs_ASDPoor_SA_effectSize.pdf"))
p

t.test(es.Good_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = 0.64854, df = 9.4989, p-value = 0.532
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2769209  0.5020282
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.07141322             -0.04114044
dotSize = 6
fontSize=28

fdr_thresh = 0.01
output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
output_res_fdr[order(output_res_fdr$pval),]
##                                    Fstat         pval      etasq         fdrq
## CortexVol                      18.768424 3.563168e-08 0.06099610 1.817216e-06
## anteromedialtemporal.lh.area   14.008061 2.095489e-06 0.10186263 5.343498e-05
## SAtotal                        12.791873 6.069546e-06 0.03840312 1.031823e-04
## superiorparietal.rh.area       10.622118 4.214024e-05 0.04068863 5.372881e-04
## posterolateraltemporal.rh.area 10.176880 6.297412e-05 0.05329274 6.423360e-04
## inferiorparietal.lh.thickness   9.067289 1.726457e-04 0.08552511 1.467488e-03
## superiorparietal.lh.area        8.890045 2.030244e-04 0.04696264 1.479178e-03
## middletemporal.rh.thickness     7.237711 9.322992e-04 0.05678581 5.943408e-03
## dorsomedialfrontal.lh.area      6.920743 1.252445e-03 0.05009390 7.097189e-03
##                                es.TD_vs_Good es.TD_vs_Poor es.Good_vs_Poor
## CortexVol                        -0.16818550    -0.6274776      -0.4168713
## anteromedialtemporal.lh.area      0.07390819    -0.7161301      -0.7587454
## SAtotal                          -0.08607947    -0.4711913      -0.3729846
## superiorparietal.rh.area         -0.23926826     0.2664748       0.5470558
## posterolateraltemporal.rh.area    0.16622170    -0.4132726      -0.6045173
## inferiorparietal.lh.thickness     0.40205640    -0.3760386      -0.7441284
## superiorparietal.lh.area         -0.33883360     0.2234276       0.5287197
## middletemporal.rh.thickness      -0.60232551    -0.3351904       0.2922828
## dorsomedialfrontal.lh.area       -0.41358040     0.1177013       0.5384547
##                                tstat.TD_vs_Good tstat.TD_vs_Poor
## CortexVol                            -0.9276362       -3.0592624
## anteromedialtemporal.lh.area          0.3724331       -3.5151991
## SAtotal                              -0.5306950       -2.1315649
## superiorparietal.rh.area             -1.5620062        1.9844711
## posterolateraltemporal.rh.area        1.1544523       -2.0474774
## inferiorparietal.lh.thickness         2.5259861       -2.1316333
## superiorparietal.lh.area             -2.0021698        1.5518600
## middletemporal.rh.thickness          -3.2211807       -1.7780643
## dorsomedialfrontal.lh.area           -2.3678548        0.5099888
##                                tstat.Good_vs_Poor pval.TD_vs_Good
## CortexVol                                2.510154     0.355355905
## anteromedialtemporal.lh.area             4.064835     0.710196116
## SAtotal                                  2.190035     0.596556966
## superiorparietal.rh.area                -2.773226     0.120794723
## posterolateraltemporal.rh.area           3.546670     0.250499638
## inferiorparietal.lh.thickness            4.082386     0.012776791
## superiorparietal.lh.area                -2.781950     0.047413637
## middletemporal.rh.thickness             -1.605353     0.001624612
## dorsomedialfrontal.lh.area              -3.412330     0.019412964
##                                pval.TD_vs_Poor pval.Good_vs_Poor
## CortexVol                         0.0026877095      1.334511e-02
## anteromedialtemporal.lh.area      0.0006040535      8.475100e-05
## SAtotal                           0.0348955186      3.037760e-02
## superiorparietal.rh.area          0.0492929914      6.408005e-03
## posterolateraltemporal.rh.area    0.0426100630      5.514425e-04
## inferiorparietal.lh.thickness     0.0349040119      7.930211e-05
## superiorparietal.lh.area          0.1231083165      6.248267e-03
## middletemporal.rh.thickness       0.0777132171      1.109600e-01
## dorsomedialfrontal.lh.area        0.6109174538      8.703711e-04
##                                fdrq.TD_vs_Good fdrq.TD_vs_Poor
## CortexVol                           0.69704428     0.022331915
## anteromedialtemporal.lh.area        0.96215984     0.007701682
## SAtotal                             0.87414672     0.148342050
## superiorparietal.rh.area            0.35564245     0.179567326
## posterolateraltemporal.rh.area      0.58419158     0.167162555
## inferiorparietal.lh.thickness       0.12381723     0.148342050
## superiorparietal.lh.area            0.19408551     0.330448639
## middletemporal.rh.thickness         0.02761841     0.247710879
## dorsomedialfrontal.lh.area          0.12381723     0.788600212
##                                fdrq.Good_vs_Poor
## CortexVol                            0.075622317
## anteromedialtemporal.lh.area         0.002161151
## SAtotal                              0.140841621
## superiorparietal.rh.area             0.046686893
## posterolateraltemporal.rh.area       0.009374523
## inferiorparietal.lh.thickness        0.002161151
## superiorparietal.lh.area             0.046686893
## middletemporal.rh.thickness          0.282947979
## dorsomedialfrontal.lh.area           0.011097232
df2plot = output_res[4:nrow(output_res),]
df2plot$region = rownames(df2plot)
new_region_labels = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal")
hemi_labels = c("LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH")
metric_labels = c("SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT")

df2plot$region = new_region_labels
df2plot$hemisphere = hemi_labels
df2plot$metric = metric_labels


colors2use = get_ggColorHue(7)

# SA F-stat
d2plot = subset(df2plot, df2plot$metric=="SA" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))

# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot, 
         display_numbers = round(mat2plot,digits=2), 
         number_color = "black", fontsize_number = 18,
         show_rownames=TRUE,
         labels_col = colnames(mat2plot),
         color = colorRampPalette(c('light blue','white','red'))(100),
         cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
                c("#785e00","#c49a00","#ffcc12"),
                c("#306800","#53b400","#76ff01"),
                c("#007459","#00c094","#0effc8"),
                c("#007b9f","#00b6eb","#39d2ff"),
                c("#6a3eff","#a58aff","#e0d7ff"),
                c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(0.6,1.3), c(0.8,1.15), c(0.75,1.35), c(0.8,1.35),
             c(0.8,1.4), c(0.8, 1.15), c(0.7,1.1))
                
for (ireg in 1:length(sig_reg)){
  region = sig_reg[ireg]
  reg_name = reg_names[ireg]
  metric_name = metric_names[ireg]
  p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
  p1 = p1 + geom_jitter(size = dotSize) +
    geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
    scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) + 
    xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) + 
    theme(text = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
  print(p1)
  
  tmp_df = d2plot[region,]
  if (tmp_df$hemisphere=="LH"){
    HEMI = "left"
  } else if (tmp_df$hemisphere=="RH"){
    HEMI = "right"
  }
  p = ggseg(.data = tmp_df, 
            atlas = chenAr, 
            mapping=aes(fill=label2fill), 
            position="stacked", 
            colour="black",
            hemisphere=HEMI) +
    scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) + 
    theme(text = element_text(size=fontSize))
  print(p)
}

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CortexVol"
reg_name = "Total Cortical Volume"
metric_name = "Volume"
ylims = c(450000,700000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
p1

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "SAtotal"
reg_name = "Total Surface Area"
metric_name = "SA"
ylims = c(90000,220000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
p1

# SA es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# SA es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# SA es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT F-stat
d2plot = subset(df2plot, df2plot$metric=="CT" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))

# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot, 
         display_numbers = round(mat2plot,digits=2), 
         number_color = "black", fontsize_number = 18,
         show_rownames=TRUE,
         labels_col = colnames(mat2plot),
         color = colorRampPalette(c('light blue','white','red'))(100),
         cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_manual(values = colors2use[1:4]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_manual(values = colors2use[5:7]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
                c("#785e00","#c49a00","#ffcc12"),
                c("#306800","#53b400","#76ff01"),
                c("#007459","#00c094","#0effc8"),
                c("#007b9f","#00b6eb","#39d2ff"),
                c("#6a3eff","#a58aff","#e0d7ff"),
                c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(-0.4,0.1), c(-0.4,0.4), c(-0.4,0.3), c(0,0.7), 
             c(-0.3,0.3), c(-0.2,0.5), c(0,0.8))
                
for (ireg in 1:length(sig_reg)){
  region = sig_reg[ireg]
  reg_name = reg_names[ireg]
  metric_name = metric_names[ireg]
  p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
  p1 = p1 + geom_jitter(size = dotSize) +
    geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
    scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) + 
    xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) + 
    theme(text = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
  print(p1)
  
  tmp_df = d2plot[region,]
  if (tmp_df$hemisphere=="LH"){
    HEMI = "left"
  } else if (tmp_df$hemisphere=="RH"){
    HEMI = "right"
  }
  p = ggseg(.data = tmp_df, 
            atlas = chenTh, 
            mapping=aes(fill=label2fill), 
            position="stacked", 
            colour="black",
            hemisphere=HEMI) +
    scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) + 
    theme(text = element_text(size=fontSize))
  print(p)
}

# CT es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

Do modeling of local differences without adjustment for total cortical volume.

# compute total SA and CT mean
asd_mri = read.csv(file.path(datapath,"asd_mri_data.csv"))
td_mri = read.csv(file.path(datapath,"td_mri_data.csv"))
lh_sa_cols = 2:35
rh_sa_cols = 106:139
lh_ct_cols = 37:70
rh_ct_cols = 141:174
asd_mri$lh_SAtotal = rowSums(asd_mri[,lh_sa_cols])
asd_mri$rh_SAtotal = rowSums(asd_mri[,rh_sa_cols])
asd_mri$SAtotal = asd_mri$lh_SAtotal + asd_mri$rh_SAtotal
asd_mri$lh_CTmean = rowMeans(asd_mri[,lh_ct_cols])
asd_mri$rh_CTmean = rowMeans(asd_mri[,rh_ct_cols])
asd_mri$CTmean = rowMeans(asd_mri[,c(lh_ct_cols,rh_ct_cols)])
td_mri$lh_SAtotal = rowSums(td_mri[,lh_sa_cols])
td_mri$rh_SAtotal = rowSums(td_mri[,rh_sa_cols])
td_mri$SAtotal = td_mri$lh_SAtotal + td_mri$rh_SAtotal
td_mri$lh_CTmean = rowMeans(td_mri[,lh_ct_cols])
td_mri$rh_CTmean = rowMeans(td_mri[,rh_ct_cols])
td_mri$CTmean = rowMeans(td_mri[,c(lh_ct_cols,rh_ct_cols)])
all_mri = rbind(td_mri,asd_mri)
all_mri$subjectId = factor(all_mri$X)

fs_data = read.csv(file.path(datapath,"FSparc_GenTemp_all.csv"))
fs_data$Group = factor(fs_data$Group)
region_names = colnames(fs_data)[5:ncol(fs_data)]
region_names = c("CortexVol","SAtotal","CTmean",region_names)

fs_data = merge(fs_data, pheno_data_mri[,c("subjectId","sex","scan_age")], by="subjectId")
fs_data = merge(fs_data, all_mri[,c("subjectId","SAtotal","CTmean")], by="subjectId")

new_fs_data = fs_data

colnames2use = c("Fstat","pval","etasq","fdrq","es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor",
                 "tstat.TD_vs_Good","tstat.TD_vs_Poor","tstat.Good_vs_Poor",
                 "pval.TD_vs_Good","pval.TD_vs_Poor","pval.Good_vs_Poor")
output_res = data.frame(matrix(nrow=length(region_names), ncol=length(colnames2use)))
rownames(output_res) = region_names
colnames(output_res) = colnames2use

for (icol in 1:length(region_names)){
  
  y_var = region_names[icol]
  
  # make formula
  form2use = as.formula(sprintf("%s ~ Group + sex + scan_age", y_var))
  full_model = model.matrix(~0+ Group + sex + scan_age, data = fs_data)
  covname2use = c("sexMale","scan_age")

  # run linear model
  mod2use = lm(formula = form2use, data = fs_data)
  
  # run an ANOVA
  res2use = anova(mod2use)
  
  # grab ANOVA results and store into output_res
  output_res[y_var, "Fstat"] = res2use["Group", "F value"]
  output_res[y_var, "pval"] = res2use["Group", "Pr(>F)"]
  
  # compute partial eta-squared as the effect size of the interaction effect
  eta_sq_res = etasq(mod2use)
  output_res[icol, "etasq"] = eta_sq_res["Group","Partial eta^2"]
  
  
    # remove covariates before effect size computation
    beta1 = mod2use$coefficients[covname2use, drop = FALSE]
    beta1[is.na(beta1)] = 0

    new_fs_data[,y_var] = as.numeric(t(fs_data[,y_var] - beta1 %*% t(full_model[,covname2use])))
    
    # compute effect sizes
    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Good"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Good"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Good"] = res$coefficients[2,"Pr(>|t|)"] 
    
    mask1 = fs_data$Group=="Poor"
    mask2 = fs_data$Group=="TD"
    output_res[y_var, "es.TD_vs_Poor"] = cohens_d(new_fs_data[mask2,y_var], new_fs_data[mask1,y_var])
    
    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.TD_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.TD_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 


    mask1 = fs_data$Group=="Good"
    mask2 = fs_data$Group=="Poor"
    output_res[y_var, "es.Good_vs_Poor"] = cohens_d(new_fs_data[mask1,y_var], new_fs_data[mask2,y_var])

    d2use = subset(fs_data, mask1 | mask2)
  mod2use = lm(formula = form2use, data = d2use)
  res = summary(mod2use)
    output_res[y_var, "tstat.Good_vs_Poor"] = res$coefficients[2,"t value"] 
    output_res[y_var, "pval.Good_vs_Poor"] = res$coefficients[2,"Pr(>|t|)"] 

}

output_res$fdrq = p.adjust(output_res$pval, method = "fdr")
output_res$fdrq.TD_vs_Good = p.adjust(output_res$pval.TD_vs_Good, method = "fdr")
output_res$fdrq.TD_vs_Poor = p.adjust(output_res$pval.TD_vs_Poor, method = "fdr")
output_res$fdrq.Good_vs_Poor = p.adjust(output_res$pval.Good_vs_Poor, method = "fdr")
output_res[order(output_res$pval),]
##                                          Fstat         pval        etasq
## CortexVol                          18.76842425 3.563168e-08 6.099610e-02
## anteromedialtemporal.lh.area       13.96198552 2.171932e-06 1.187908e-01
## SAtotal                            12.79187344 6.069546e-06 3.840312e-02
## superiorparietal.rh.area            9.97848940 7.518206e-05 6.112751e-02
## posterolateraltemporal.rh.area      9.79945558 8.843829e-05 7.466123e-02
## inferiorparietal.lh.thickness       9.11412690 1.650865e-04 8.726838e-02
## superiorparietal.lh.area            8.39947078 3.178795e-04 6.424037e-02
## middletemporal.rh.thickness         7.25613913 9.152796e-04 5.855042e-02
## dorsomedialfrontal.lh.area          6.95235733 1.214613e-03 4.948907e-02
## perisylvian.lh.thickness            5.52500358 4.641913e-03 2.703457e-02
## inferiorparietal.rh.thickness       5.13340268 6.727702e-03 4.630736e-02
## middletemporal.lh.thickness         4.88755620 8.498974e-03 3.322598e-02
## orbitofrontal.rh.area               4.73445787 9.833290e-03 8.651774e-02
## dorsomedialfrontal.rh.area          4.67781875 1.037894e-02 4.928688e-02
## medialprefrontal.rh.thickness       4.26035165 1.546793e-02 1.189294e-02
## medialprefrontal.lh.thickness       4.12360606 1.763387e-02 1.682454e-02
## inferiorparietal.rh.area            4.05447261 1.884306e-02 4.356279e-02
## orbitofrontal.lh.area               3.97768109 2.028472e-02 7.394201e-02
## occipital.lh.area                   3.54678869 3.071012e-02 2.442052e-02
## posterolateraltemporal.lh.area      3.47096883 3.304108e-02 2.333083e-02
## dorsolateralprefrotal.lh.thickness  3.38761197 3.581075e-02 1.222666e-02
## temporalpole.lh.thickness           3.10692393 4.698360e-02 2.235669e-02
## ventromedialoccipital.lh.thickness  2.84207875 6.074741e-02 4.976358e-02
## superiorparietal.rh.thickness       2.80424833 6.302173e-02 2.572353e-02
## anteromedialtemporal.rh.area        2.70365990 6.949654e-02 8.384856e-03
## medialtemporal.lh.thickness         2.57867720 7.848677e-02 2.562891e-02
## CTmean                              2.52041363 8.307086e-02 1.158927e-03
## inferiorparietal.lh.area            2.51632609 8.340244e-02 2.073572e-02
## dorsolateralprefrotal.rh.thickness  2.43567482 9.022542e-02 6.603310e-03
## temporalpole.rh.thickness           2.25922675 1.071864e-01 9.199317e-03
## perisylvian.rh.thickness            2.16026648 1.180743e-01 1.186715e-02
## motor_premotor_SMA.rh.thickness     2.10997052 1.240297e-01 4.047085e-02
## superiortemporal.rh.area            2.06577098 1.295134e-01 5.537640e-03
## ventralfrontal.rh.thickness         1.98772991 1.398011e-01 1.206705e-02
## motor_premotor_SMA.lh.thickness     1.61382458 2.017999e-01 5.746284e-02
## dorsolateralprefrontal.rh.area      1.38860475 2.519036e-01 1.261047e-02
## occipital.lh.thickness              1.28537276 2.789040e-01 1.422125e-02
## occipital.rh.thickness              1.26440340 2.847364e-01 1.084307e-02
## precuneus.rh.area                   1.25553820 2.872390e-01 1.159653e-02
## motorpremotor.rh.area               1.09409666 3.369095e-01 2.026216e-02
## ventromedialoccipital.rh.thickness  1.08302969 3.406168e-01 1.771466e-02
## parsopercularis.lh.area             1.02483692 3.607994e-01 7.390424e-03
## ventralfrontal.lh.thickness         0.83696514 4.345882e-01 6.276432e-03
## dorsolateralprefrontal.lh.area      0.53335542 5.874954e-01 8.638171e-03
## precuneus.lh.area                   0.48186041 6.183743e-01 8.298034e-03
## medialtemporal.rh.thickness         0.27931431 7.566073e-01 3.548265e-03
## motorpremotor.lh.area               0.26497128 7.675066e-01 1.317870e-03
## occipital.rh.area                   0.20422404 8.154556e-01 2.499252e-03
## parsopercularis.rh.area             0.17081341 8.431061e-01 1.638071e-04
## superiortemporal.lh.area            0.07917468 9.239085e-01 2.428987e-03
## superiorparietal.lh.thickness       0.02132311 9.789049e-01 2.247971e-05
##                                            fdrq es.TD_vs_Good es.TD_vs_Poor
## CortexVol                          1.817216e-06  -0.168185504  -0.627477606
## anteromedialtemporal.lh.area       5.538426e-05   0.058564872  -0.779840965
## SAtotal                            1.031823e-04  -0.086079474  -0.471191290
## superiorparietal.rh.area           9.020706e-04  -0.192365758   0.422061004
## posterolateraltemporal.rh.area     9.020706e-04   0.132202091  -0.530611464
## inferiorparietal.lh.thickness      1.403235e-03   0.403204612  -0.371601402
## superiorparietal.lh.area           2.315979e-03  -0.284174733   0.382993616
## middletemporal.rh.thickness        5.834907e-03  -0.607526883  -0.368109847
## dorsomedialfrontal.lh.area         6.882808e-03  -0.417777838   0.101394689
## perisylvian.lh.thickness           2.367376e-02  -0.417755510  -0.141804387
## inferiorparietal.rh.thickness      3.119207e-02   0.436965414  -0.080961775
## middletemporal.lh.thickness        3.612064e-02  -0.111703511  -0.484911933
## orbitofrontal.rh.area              3.780900e-02   0.693881222   0.663368868
## dorsomedialfrontal.rh.area         3.780900e-02  -0.017052068   0.536607213
## medialprefrontal.rh.thickness      5.259096e-02   0.108471803   0.272612686
## medialprefrontal.lh.thickness      5.620796e-02   0.034413715   0.303600784
## inferiorparietal.rh.area           5.652917e-02  -0.465146831   0.088250958
## orbitofrontal.lh.area              5.747337e-02   0.302307914   0.684921768
## occipital.lh.area                  8.243242e-02   0.405773203   0.206264627
## posterolateraltemporal.lh.area     8.425477e-02   0.040680967  -0.313873302
## dorsolateralprefrotal.lh.thickness 8.696896e-02   0.100397423   0.287923700
## temporalpole.lh.thickness          1.089165e-01  -0.036104455  -0.356907842
## ventromedialoccipital.lh.thickness 1.339212e-01  -0.426319170  -0.598130176
## superiorparietal.rh.thickness      1.339212e-01   0.287207325  -0.083702106
## anteromedialtemporal.rh.area       1.417729e-01  -0.068797641  -0.236087202
## medialtemporal.lh.thickness        1.519116e-01   0.086741514   0.393647416
## CTmean                             1.519116e-01   0.022296408   0.086644973
## inferiorparietal.lh.area           1.519116e-01  -0.264343885  -0.361967769
## dorsolateralprefrotal.rh.thickness 1.586723e-01   0.158394449  -0.019007650
## temporalpole.rh.thickness          1.822168e-01  -0.155643965  -0.252523004
## perisylvian.rh.thickness           1.942512e-01  -0.112868711  -0.270616154
## motor_premotor_SMA.rh.thickness    1.976724e-01   0.283654397   0.508152496
## superiortemporal.rh.area           2.001571e-01   0.181205381   0.133269647
## ventralfrontal.rh.thickness        2.097016e-01  -0.245560738  -0.034176376
## motor_premotor_SMA.lh.thickness    2.940513e-01   0.532800679   0.594555794
## dorsolateralprefrontal.rh.area     3.568634e-01   0.141028172  -0.131031873
## occipital.lh.thickness             3.756202e-01  -0.101511168   0.187461609
## occipital.rh.thickness             3.756202e-01   0.003233587   0.234784901
## precuneus.rh.area                  3.756202e-01  -0.116016165   0.142181444
## motorpremotor.rh.area              4.236940e-01  -0.219240608  -0.350517881
## ventromedialoccipital.rh.thickness 4.236940e-01  -0.309481377  -0.275857079
## parsopercularis.lh.area            4.381135e-01   0.038454375  -0.187192775
## ventralfrontal.lh.thickness        5.154418e-01  -0.201740570  -0.035264557
## dorsolateralprefrontal.lh.area     6.809605e-01   0.146391233   0.230406377
## precuneus.lh.area                  7.008242e-01  -0.055471505   0.161218164
## medialtemporal.rh.thickness        8.328263e-01   0.054942375   0.168683420
## motorpremotor.lh.area              8.328263e-01  -0.023422791   0.062490351
## occipital.rh.area                  8.664216e-01   0.056561428   0.118308443
## parsopercularis.rh.area            8.775186e-01  -0.027403626  -0.027013701
## superiortemporal.lh.area           9.423867e-01  -0.120637855  -0.104382282
## superiorparietal.lh.thickness      9.789049e-01   0.011660023   0.008979405
##                                    es.Good_vs_Poor tstat.TD_vs_Good
## CortexVol                             -0.416871265      -0.92763624
## anteromedialtemporal.lh.area          -0.794803988       0.16733131
## SAtotal                               -0.372984599      -0.53069497
## superiorparietal.rh.area               0.627981093      -1.28721394
## posterolateraltemporal.rh.area        -0.684868006       0.96665266
## inferiorparietal.lh.thickness         -0.741115244       2.47252195
## superiorparietal.lh.area               0.617690683      -1.58921178
## middletemporal.rh.thickness            0.267713645      -3.33001733
## dorsomedialfrontal.lh.area             0.525846171      -2.35106489
## perisylvian.lh.thickness               0.272324607      -2.37407644
## inferiorparietal.rh.thickness         -0.529004839       2.51701196
## middletemporal.lh.thickness           -0.320912774      -0.64723145
## orbitofrontal.rh.area                 -0.043209574       3.72905434
## dorsomedialfrontal.rh.area             0.466219000      -0.04161523
## medialprefrontal.rh.thickness          0.161735686       0.88959051
## medialprefrontal.lh.thickness          0.247253521       0.14898278
## inferiorparietal.rh.area               0.446921552      -2.83757132
## orbitofrontal.lh.area                  0.418212137       1.63638457
## occipital.lh.area                     -0.188622586       2.13738216
## posterolateraltemporal.lh.area        -0.360688505       0.09481595
## dorsolateralprefrotal.lh.thickness     0.173019231       0.46757198
## temporalpole.lh.thickness             -0.300282103      -0.24817897
## ventromedialoccipital.lh.thickness    -0.085302255      -2.38119144
## superiorparietal.rh.thickness         -0.403949755       1.77629274
## anteromedialtemporal.rh.area          -0.148614009      -0.53739563
## medialtemporal.lh.thickness            0.303841318       0.45893648
## CTmean                                 0.059512817       0.06023806
## inferiorparietal.lh.area              -0.116200804      -1.32874444
## dorsolateralprefrotal.rh.thickness    -0.196139656       0.47428576
## temporalpole.rh.thickness             -0.080405337      -0.92866917
## perisylvian.rh.thickness              -0.163250435      -0.71061831
## motor_premotor_SMA.rh.thickness        0.238945563       1.65793388
## superiortemporal.rh.area              -0.045635530       1.32733836
## ventralfrontal.rh.thickness            0.218413802      -1.59784005
## motor_premotor_SMA.lh.thickness        0.034411163       3.10024557
## dorsolateralprefrontal.rh.area        -0.297885547       0.89205477
## occipital.lh.thickness                 0.305012767      -0.66781162
## occipital.rh.thickness                 0.219594316       0.02112775
## precuneus.rh.area                      0.295072659      -0.64623092
## motorpremotor.rh.area                 -0.139743512      -1.13212189
## ventromedialoccipital.rh.thickness     0.068451005      -1.66827098
## parsopercularis.lh.area               -0.195904940       0.12501934
## ventralfrontal.lh.thickness            0.142797772      -1.35421352
## dorsolateralprefrontal.lh.area         0.088096781       0.92201234
## precuneus.lh.area                      0.222805118      -0.19224401
## medialtemporal.rh.thickness            0.086489168       0.38406302
## motorpremotor.lh.area                  0.086177916      -0.24000071
## occipital.rh.area                      0.071541806       0.13971882
## parsopercularis.rh.area                0.003238283      -0.09907955
## superiortemporal.lh.area               0.016663201      -0.21163657
## superiorparietal.lh.thickness         -0.003378940       0.27808249
##                                    tstat.TD_vs_Poor tstat.Good_vs_Poor
## CortexVol                              -3.059262385         2.51015442
## anteromedialtemporal.lh.area           -3.457718886         4.44548917
## SAtotal                                -2.131564929         2.19003479
## superiorparietal.rh.area                2.459248096        -3.58518048
## posterolateraltemporal.rh.area         -2.687567084         4.10526246
## inferiorparietal.lh.thickness          -2.045168291         4.13088855
## superiorparietal.lh.area                2.027310188        -3.41805335
## middletemporal.rh.thickness            -1.747812746        -1.47491113
## dorsomedialfrontal.lh.area              0.493902033        -3.20521339
## perisylvian.lh.thickness               -0.804089909        -1.53042696
## inferiorparietal.rh.thickness          -1.006601493         2.65923159
## middletemporal.lh.thickness            -2.249001377         1.86010176
## orbitofrontal.rh.area                   4.159885029         0.32847869
## dorsomedialfrontal.rh.area              3.102459847        -2.54332578
## medialprefrontal.rh.thickness           1.200505533        -0.82509176
## medialprefrontal.lh.thickness           1.480004889        -1.40025248
## inferiorparietal.rh.area                0.661789105        -2.43272506
## orbitofrontal.lh.area                   3.660090310        -2.40201886
## occipital.lh.area                       0.791523855         0.99772299
## posterolateraltemporal.lh.area         -1.116709014         2.19083651
## dorsolateralprefrotal.lh.thickness      1.511150424        -0.92199454
## temporalpole.lh.thickness              -1.870933210         1.82402736
## ventromedialoccipital.lh.thickness     -3.457098374         0.21658321
## superiorparietal.rh.thickness          -0.406187547         2.20945034
## anteromedialtemporal.rh.area           -0.952204003         0.83690241
## medialtemporal.lh.thickness             2.192865606        -1.44986930
## CTmean                                  0.003809356        -0.36930388
## inferiorparietal.lh.area               -2.167118278         0.72881076
## dorsolateralprefrotal.rh.thickness      0.019698750         0.99601950
## temporalpole.rh.thickness              -1.230129724         0.75667370
## perisylvian.rh.thickness               -1.517818628         0.77574387
## motor_premotor_SMA.rh.thickness         2.832075777        -1.28000273
## superiortemporal.rh.area                0.112341163         0.34644699
## ventralfrontal.rh.thickness            -0.549063137        -1.31651197
## motor_premotor_SMA.lh.thickness         3.517773482        -0.12326995
## dorsolateralprefrontal.rh.area         -0.610137269         1.68378954
## occipital.lh.thickness                  0.546421993        -1.93570525
## occipital.rh.thickness                  1.235174967        -1.27724201
## precuneus.rh.area                       1.048263005        -1.84337987
## motorpremotor.rh.area                  -2.128265711         0.67938285
## ventromedialoccipital.rh.thickness     -0.804993424        -0.33280224
## parsopercularis.lh.area                -1.396669225         0.99693457
## ventralfrontal.lh.thickness            -0.236705025        -0.77527381
## dorsolateralprefrontal.lh.area          0.943125104        -0.51787823
## precuneus.lh.area                       1.148868563        -1.27347487
## medialtemporal.rh.thickness             1.351096705        -0.04991803
## motorpremotor.lh.area                   0.669118895        -0.41733396
## occipital.rh.area                       0.260126462        -0.55310296
## parsopercularis.rh.area                -0.598475549        -0.13375612
## superiortemporal.lh.area               -1.371031861        -0.10594187
## superiorparietal.lh.thickness           0.180689829        -0.04619990
##                                    pval.TD_vs_Good pval.TD_vs_Poor
## CortexVol                             0.3553559051    0.0026877095
## anteromedialtemporal.lh.area          0.8673755167    0.0007336348
## SAtotal                               0.5965569661    0.0348955186
## superiorparietal.rh.area              0.2003607435    0.0152169867
## posterolateraltemporal.rh.area        0.3355548864    0.0081243665
## inferiorparietal.lh.thickness         0.0147388815    0.0428253247
## superiorparietal.lh.area              0.1144990012    0.0446439975
## middletemporal.rh.thickness           0.0011370894    0.0828217214
## dorsomedialfrontal.lh.area            0.0202592759    0.6221962228
## perisylvian.lh.thickness              0.0190909801    0.4227907979
## inferiorparietal.rh.thickness         0.0130803984    0.3159679753
## middletemporal.lh.thickness           0.5186503403    0.0261692628
## orbitofrontal.rh.area                 0.0002883625    0.0000569994
## dorsomedialfrontal.rh.area            0.9668707826    0.0023477810
## medialprefrontal.rh.thickness         0.3753677892    0.2320930412
## medialprefrontal.lh.thickness         0.8818035363    0.1412547090
## inferiorparietal.rh.area              0.0052934622    0.5092600271
## orbitofrontal.lh.area                 0.1042349448    0.0003634346
## occipital.lh.area                     0.0344834500    0.4300583771
## posterolateraltemporal.lh.area        0.9246104541    0.2661473395
## dorsolateralprefrotal.lh.thickness    0.6408920214    0.1331404116
## temporalpole.lh.thickness             0.8043967759    0.0635682953
## ventromedialoccipital.lh.thickness    0.0187420175    0.0007351844
## superiorparietal.rh.thickness         0.0780791851    0.6852625412
## anteromedialtemporal.rh.area          0.5919347149    0.3427338726
## medialtemporal.lh.thickness           0.6470642744    0.0300690708
## CTmean                                0.9520607209    0.9969663319
## inferiorparietal.lh.area              0.1863147768    0.0320206760
## dorsolateralprefrotal.rh.thickness    0.6361105984    0.9843134327
## temporalpole.rh.thickness             0.3548222583    0.2208363467
## perisylvian.rh.thickness              0.4786237488    0.1314515201
## motor_premotor_SMA.rh.thickness       0.0997988784    0.0053501067
## superiortemporal.rh.area              0.1867779598    0.9107235653
## ventralfrontal.rh.thickness           0.1125637166    0.5838897264
## motor_premotor_SMA.lh.thickness       0.0023821153    0.0005973643
## dorsolateralprefrontal.rh.area        0.3740506397    0.5428198559
## occipital.lh.thickness                0.5054656127    0.5856983221
## occipital.rh.thickness                0.9831769052    0.2189593711
## precuneus.rh.area                     0.5192958825    0.2964330144
## motorpremotor.rh.area                 0.2597162243    0.0351732624
## ventromedialoccipital.rh.thickness    0.0977256551    0.4222710686
## parsopercularis.lh.area               0.9007059164    0.1648569587
## ventralfrontal.lh.thickness           0.1780729303    0.8132525724
## dorsolateralprefrontal.lh.area        0.3582703860    0.3473398905
## precuneus.lh.area                     0.8478581339    0.2526877319
## medialtemporal.rh.thickness           0.7015740123    0.1789757279
## motorpremotor.lh.area                 0.8107166551    0.5045882392
## occipital.rh.area                     0.8891034877    0.7951715449
## parsopercularis.rh.area               0.9212313305    0.5505482448
## superiortemporal.lh.area              0.8327298285    0.1726915381
## superiorparietal.lh.thickness         0.7814012556    0.8568881688
##                                    pval.Good_vs_Poor fdrq.TD_vs_Good
## CortexVol                               1.334511e-02      0.70902805
## anteromedialtemporal.lh.area            1.911862e-05      0.98239861
## SAtotal                                 3.037760e-02      0.91667439
## superiorparietal.rh.area                4.814179e-04      0.51091990
## posterolateraltemporal.rh.area          7.238701e-05      0.70902805
## inferiorparietal.lh.thickness           6.564083e-05      0.11480256
## superiorparietal.lh.area                8.518760e-04      0.36496557
## middletemporal.rh.thickness             1.427507e-01      0.02899578
## dorsomedialfrontal.lh.area              1.713282e-03      0.11480256
## perisylvian.lh.thickness                1.284379e-01      0.11480256
## inferiorparietal.rh.thickness           8.857225e-03      0.11480256
## middletemporal.lh.thickness             6.522152e-02      0.85432548
## orbitofrontal.rh.area                   7.430994e-01      0.01470649
## dorsomedialfrontal.rh.area              1.219960e-02      0.98317691
## medialprefrontal.rh.thickness           4.108913e-01      0.70902805
## medialprefrontal.lh.thickness           1.639150e-01      0.98239861
## inferiorparietal.rh.area                1.640132e-02      0.06749164
## orbitofrontal.lh.area                   1.777626e-02      0.36496557
## occipital.lh.area                       3.203412e-01      0.17586560
## posterolateraltemporal.lh.area          3.031810e-02      0.98239861
## dorsolateralprefrotal.lh.thickness      3.583076e-01      0.91667439
## temporalpole.lh.thickness               7.053658e-02      0.98239861
## ventromedialoccipital.lh.thickness      8.288862e-01      0.11480256
## superiorparietal.rh.thickness           2.896481e-02      0.36200349
## anteromedialtemporal.rh.area            4.042446e-01      0.91667439
## medialtemporal.lh.thickness             1.495994e-01      0.91667439
## CTmean                                  7.125261e-01      0.98317691
## inferiorparietal.lh.area                4.674807e-01      0.50135137
## dorsolateralprefrotal.rh.thickness      3.211649e-01      0.91667439
## temporalpole.rh.thickness               4.506692e-01      0.70902805
## perisylvian.rh.thickness                4.393654e-01      0.85432548
## motor_premotor_SMA.rh.thickness         2.029146e-01      0.36496557
## superiortemporal.rh.area                7.295893e-01      0.50135137
## ventralfrontal.rh.thickness             1.904111e-01      0.36496557
## motor_premotor_SMA.lh.thickness         9.020913e-01      0.04049596
## dorsolateralprefrontal.rh.area          9.471804e-02      0.70902805
## occipital.lh.thickness                  5.516167e-02      0.85432548
## occipital.rh.thickness                  2.038841e-01      0.98317691
## precuneus.rh.area                       6.764219e-02      0.85432548
## motorpremotor.rh.area                   4.981514e-01      0.63073940
## ventromedialoccipital.rh.thickness      7.398411e-01      0.36496557
## parsopercularis.lh.area                 3.207222e-01      0.98239861
## ventralfrontal.lh.thickness             4.396421e-01      0.50135137
## dorsolateralprefrontal.lh.area          6.054586e-01      0.70902805
## precuneus.lh.area                       2.052126e-01      0.98239861
## medialtemporal.rh.thickness             9.602674e-01      0.96703445
## motorpremotor.lh.area                   6.771501e-01      0.98239861
## occipital.rh.area                       5.811808e-01      0.98239861
## parsopercularis.rh.area                 8.938105e-01      0.98239861
## superiortemporal.lh.area                9.157983e-01      0.98239861
## superiorparietal.lh.thickness           9.632247e-01      0.98239861
##                                    fdrq.TD_vs_Poor fdrq.Good_vs_Poor
## CortexVol                              0.019581884      0.0756223169
## anteromedialtemporal.lh.area           0.007498881      0.0009750497
## SAtotal                                0.119589092      0.1106612739
## superiorparietal.rh.area               0.077606632      0.0061380785
## posterolateraltemporal.rh.area         0.046038077      0.0012305792
## inferiorparietal.lh.thickness          0.133931993      0.0012305792
## superiorparietal.lh.area               0.133931993      0.0086891350
## middletemporal.rh.thickness            0.222310936      0.3466803416
## dorsomedialfrontal.lh.area             0.721181985      0.0145628942
## perisylvian.lh.thickness               0.592783168      0.3275166530
## inferiorparietal.rh.thickness          0.503573961      0.0645312131
## middletemporal.lh.thickness            0.119589092      0.1998536554
## orbitofrontal.rh.area                  0.002906969      0.8421793049
## dorsomedialfrontal.rh.area             0.019581884      0.0756223169
## medialprefrontal.rh.thickness          0.422740896      0.6350138435
## medialprefrontal.lh.thickness          0.327454098      0.3634636481
## inferiorparietal.rh.area               0.665955420      0.0824171823
## orbitofrontal.lh.area                  0.007498881      0.0824171823
## occipital.lh.area                      0.592783168      0.5459802676
## posterolateraltemporal.lh.area         0.452450477      0.1106612739
## dorsolateralprefrotal.lh.thickness     0.323341000      0.5894737178
## temporalpole.lh.thickness              0.180110170      0.1998536554
## ventromedialoccipital.lh.thickness     0.007498881      0.9189825748
## superiorparietal.rh.thickness          0.776630880      0.1106612739
## anteromedialtemporal.rh.area           0.521009836      0.6350138435
## medialtemporal.lh.thickness            0.119589092      0.3467986259
## CTmean                                 0.996966332      0.8421793049
## inferiorparietal.lh.area               0.119589092      0.6443652932
## dorsolateralprefrotal.rh.thickness     0.996966332      0.5459802676
## temporalpole.rh.thickness              0.417135322      0.6384480199
## perisylvian.rh.thickness               0.323341000      0.6384480199
## motor_premotor_SMA.rh.thickness        0.034106930      0.3876238194
## superiortemporal.rh.area               0.947895956      0.8421793049
## ventralfrontal.rh.thickness            0.694665452      0.3876238194
## motor_premotor_SMA.lh.thickness        0.007498881      0.9531778237
## dorsolateralprefrontal.rh.area         0.684828305      0.2542431500
## occipital.lh.thickness                 0.694665452      0.1875496803
## occipital.rh.thickness                 0.417135322      0.3876238194
## precuneus.rh.area                      0.487680120      0.1998536554
## motorpremotor.rh.area                  0.119589092      0.6685716490
## ventromedialoccipital.rh.thickness     0.592783168      0.8421793049
## parsopercularis.lh.area                0.365110485      0.5459802676
## ventralfrontal.lh.thickness            0.882465557      0.6384480199
## dorsolateralprefrontal.lh.area         0.521009836      0.7719597195
## precuneus.lh.area                      0.444381873      0.3876238194
## medialtemporal.rh.thickness            0.365110485      0.9632246542
## motorpremotor.lh.area                  0.665955420      0.8421793049
## occipital.rh.area                      0.881603235      0.7600056174
## parsopercularis.rh.area                0.684828305      0.9531778237
## superiortemporal.lh.area               0.365110485      0.9531778237
## superiorparietal.lh.thickness          0.910443679      0.9632246542
regions2use = colnames(fs_data)[5:(ncol(fs_data)-4)]
full_regnames = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor","occipital","posterolateral temporal","superior parietal",
                  "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                  "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal",
                  "motor-premotor-supplementary motor area","ventral frontal","ventromedial occipital","inferior parietal",
                  "middle temporal","perisylvian","dorsolateral prefrontal","occipital",
                  "temporal pole","superior parietal","medial prefrontal","medial temporal")

sa_ordering = c(1,12,7,10,
                4,6,9,3,
                8,11,2,5)
ct_ordering = c(1,7,6,3,
                10,4,11,5,
                8,2,12,9)
gensim_order = c(sa_ordering,sa_ordering, ct_ordering,ct_ordering)

df4plot = output_res[regions2use,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")]
df4plot$phenotype = "CT"
df4plot$phenotype[1:24] = "SA"
df4plot$phenotype = factor(df4plot$phenotype) 

df4plot$hemisphere = "RH"
df4plot$hemisphere[1:12] = "LH"
df4plot$hemisphere[25:36] = "LH"
df4plot$hemisphere = factor(df4plot$hemisphere) 

df4plot$full_region_names = factor(full_regnames)
df4plot$gensim_gradient = gensim_order

ct_df = subset(df4plot, df4plot$phenotype=="CT")
sa_df = subset(df4plot, df4plot$phenotype=="SA")

ct_df$cluster = "dorsal"
ct_df$cluster[c(6,7,8,9,10,11,12)] = "ventral"
ct_df$cluster = factor(ct_df$cluster)

sa_df$cluster = "anterior"
sa_df$cluster[c(6,7,8,9,10,11,12)] = "posterior"
sa_df$cluster = factor(sa_df$cluster)

dotSize = 10
fontSize = 20

p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Good
## t = -1.661, df = 22, p-value = 0.1109
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.64971696  0.08041565
## sample estimates:
##        cor 
## -0.3338096
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Good ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 0.50429, df = 19.576, p-value = 0.6197
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1627292  0.2663074
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           0.005766454          -0.046022637
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.TD_vs_Poor
## t = -0.26326, df = 22, p-value = 0.7948
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4492788  0.3553922
## sample estimates:
##         cor 
## -0.05603835
p = ggplot(data = ct_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = -1.2347, df = 13.777, p-value = 0.2376
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4439168  0.1198526
## sample estimates:
##  mean in group dorsal mean in group ventral 
##           -0.06447487            0.09755724
p = ggplot(data = ct_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("CT ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(ct_df$gensim_gradient, ct_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  ct_df$gensim_gradient and ct_df$es.Good_vs_Poor
## t = 1.3098, df = 22, p-value = 0.2038
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1507977  0.6065479
## sample estimates:
##       cor 
## 0.2689602
p = ggplot(data = ct_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Dorsal-Ventral Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.Good_vs_Poor ~ cluster, data = ct_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = -1.8651, df = 14.85, p-value = 0.08205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.44196835  0.02964677
## sample estimates:
##  mean in group dorsal mean in group ventral 
##            -0.0636193             0.1425415
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Good, fill = es.TD_vs_Good)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Good)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Good
## t = -0.5984, df = 22, p-value = 0.5557
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5042089  0.2917377
## sample estimates:
##        cor 
## -0.1265525
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Good, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Good") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Good ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Good by cluster
## t = 1.197, df = 15.567, p-value = 0.2492
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09268998  0.33188479
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.03176589             -0.08783151
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.TD_vs_Poor, fill = es.TD_vs_Poor)) +
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.TD_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.TD_vs_Poor
## t = -0.13074, df = 22, p-value = 0.8972
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4264660  0.3798029
## sample estimates:
##         cor 
## -0.02786268
p = ggplot(data = sa_df, aes(x = cluster, y = es.TD_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters TD vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.TD_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.TD_vs_Poor by cluster
## t = 1.5331, df = 11.265, p-value = 0.1529
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1049842  0.5914930
## sample estimates:
##  mean in group anterior mean in group posterior 
##               0.1089166              -0.1343378
p = ggplot(data = sa_df, aes(x = gensim_gradient, y = es.Good_vs_Poor, fill = es.Good_vs_Poor)) + 
  geom_point(size=dotSize, colour="black", pch=21) + 
  geom_smooth(method=lm, colour="black") + xlim(0,12) +
  scale_x_continuous(breaks=seq(1,12,1)) + 
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100)) +
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Brain Region") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("SA ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

cor.test(sa_df$gensim_gradient, sa_df$es.Good_vs_Poor)
## 
##  Pearson's product-moment correlation
## 
## data:  sa_df$gensim_gradient and sa_df$es.Good_vs_Poor
## t = 0.32019, df = 22, p-value = 0.7518
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3447631  0.4588948
## sample estimates:
##       cor 
## 0.0681055
p = ggplot(data = sa_df, aes(x = cluster, y = es.Good_vs_Poor, fill=cluster))
p = p + geom_jitter(size=dotSize, colour="black", pch=21) + 
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + 
  # scale_fill_manual(values = c("#1e90ff","#ff8d1e")) + 
  scale_fill_manual(values = c("blue","red")) + 
  guides(colour=FALSE, fill=FALSE) +  
  xlab("Clusters") + 
  ylab("Effect Size (Cohen's d)") + 
  ggtitle("Anterior-Posterior Clusters ASD Good vs ASD Poor") +
  theme(text = element_text(size=fontSize), 
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
p

t.test(es.Good_vs_Poor ~ cluster, data = sa_df)
## 
##  Welch Two Sample t-test
## 
## data:  es.Good_vs_Poor by cluster
## t = 0.5723, df = 10.21, p-value = 0.5795
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2949289  0.4995674
## sample estimates:
##  mean in group anterior mean in group posterior 
##              0.06610517             -0.03621407
dotSize = 6
fontSize=28

fdr_thresh = 0.01
output_res_fdr = output_res[output_res$fdrq<=fdr_thresh,]
output_res_fdr[order(output_res_fdr$pval),]
##                                    Fstat         pval      etasq         fdrq
## CortexVol                      18.768424 3.563168e-08 0.06099610 1.817216e-06
## anteromedialtemporal.lh.area   13.961986 2.171932e-06 0.11879079 5.538426e-05
## SAtotal                        12.791873 6.069546e-06 0.03840312 1.031823e-04
## superiorparietal.rh.area        9.978489 7.518206e-05 0.06112751 9.020706e-04
## posterolateraltemporal.rh.area  9.799456 8.843829e-05 0.07466123 9.020706e-04
## inferiorparietal.lh.thickness   9.114127 1.650865e-04 0.08726838 1.403235e-03
## superiorparietal.lh.area        8.399471 3.178795e-04 0.06424037 2.315979e-03
## middletemporal.rh.thickness     7.256139 9.152796e-04 0.05855042 5.834907e-03
## dorsomedialfrontal.lh.area      6.952357 1.214613e-03 0.04948907 6.882808e-03
##                                es.TD_vs_Good es.TD_vs_Poor es.Good_vs_Poor
## CortexVol                        -0.16818550    -0.6274776      -0.4168713
## anteromedialtemporal.lh.area      0.05856487    -0.7798410      -0.7948040
## SAtotal                          -0.08607947    -0.4711913      -0.3729846
## superiorparietal.rh.area         -0.19236576     0.4220610       0.6279811
## posterolateraltemporal.rh.area    0.13220209    -0.5306115      -0.6848680
## inferiorparietal.lh.thickness     0.40320461    -0.3716014      -0.7411152
## superiorparietal.lh.area         -0.28417473     0.3829936       0.6176907
## middletemporal.rh.thickness      -0.60752688    -0.3681098       0.2677136
## dorsomedialfrontal.lh.area       -0.41777784     0.1013947       0.5258462
##                                tstat.TD_vs_Good tstat.TD_vs_Poor
## CortexVol                            -0.9276362        -3.059262
## anteromedialtemporal.lh.area          0.1673313        -3.457719
## SAtotal                              -0.5306950        -2.131565
## superiorparietal.rh.area             -1.2872139         2.459248
## posterolateraltemporal.rh.area        0.9666527        -2.687567
## inferiorparietal.lh.thickness         2.4725220        -2.045168
## superiorparietal.lh.area             -1.5892118         2.027310
## middletemporal.rh.thickness          -3.3300173        -1.747813
## dorsomedialfrontal.lh.area           -2.3510649         0.493902
##                                tstat.Good_vs_Poor pval.TD_vs_Good
## CortexVol                                2.510154     0.355355905
## anteromedialtemporal.lh.area             4.445489     0.867375517
## SAtotal                                  2.190035     0.596556966
## superiorparietal.rh.area                -3.585180     0.200360743
## posterolateraltemporal.rh.area           4.105262     0.335554886
## inferiorparietal.lh.thickness            4.130889     0.014738882
## superiorparietal.lh.area                -3.418053     0.114499001
## middletemporal.rh.thickness             -1.474911     0.001137089
## dorsomedialfrontal.lh.area              -3.205213     0.020259276
##                                pval.TD_vs_Poor pval.Good_vs_Poor
## CortexVol                         0.0026877095      1.334511e-02
## anteromedialtemporal.lh.area      0.0007336348      1.911862e-05
## SAtotal                           0.0348955186      3.037760e-02
## superiorparietal.rh.area          0.0152169867      4.814179e-04
## posterolateraltemporal.rh.area    0.0081243665      7.238701e-05
## inferiorparietal.lh.thickness     0.0428253247      6.564083e-05
## superiorparietal.lh.area          0.0446439975      8.518760e-04
## middletemporal.rh.thickness       0.0828217214      1.427507e-01
## dorsomedialfrontal.lh.area        0.6221962228      1.713282e-03
##                                fdrq.TD_vs_Good fdrq.TD_vs_Poor
## CortexVol                           0.70902805     0.019581884
## anteromedialtemporal.lh.area        0.98239861     0.007498881
## SAtotal                             0.91667439     0.119589092
## superiorparietal.rh.area            0.51091990     0.077606632
## posterolateraltemporal.rh.area      0.70902805     0.046038077
## inferiorparietal.lh.thickness       0.11480256     0.133931993
## superiorparietal.lh.area            0.36496557     0.133931993
## middletemporal.rh.thickness         0.02899578     0.222310936
## dorsomedialfrontal.lh.area          0.11480256     0.721181985
##                                fdrq.Good_vs_Poor
## CortexVol                           0.0756223169
## anteromedialtemporal.lh.area        0.0009750497
## SAtotal                             0.1106612739
## superiorparietal.rh.area            0.0061380785
## posterolateraltemporal.rh.area      0.0012305792
## inferiorparietal.lh.thickness       0.0012305792
## superiorparietal.lh.area            0.0086891350
## middletemporal.rh.thickness         0.3466803416
## dorsomedialfrontal.lh.area          0.0145628942
df2plot = output_res[4:nrow(output_res),]
df2plot$region = rownames(df2plot)
new_region_labels = c("motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor","occipital","posterolateral temporal","superior parietal",
                      "orbitofrontal","superior temporal","inferior parietal","dorsomedial frontal",
                      "anteromedial temporal","precuneus","dorsolateral prefrontal","pars opercularis,subcentral",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal",
                      "motor-premotor-supplementary motor area","ventral frontal",
                      "ventromedial occipital","inferior parietal","middle temporal","perisylvian",
                      "dorsolateral prefrontal","occipital","temporal pole","superior parietal",
                      "medial prefrontal","medial temporal")
hemi_labels = c("LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "LH","LH","LH","LH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH",
                      "RH","RH","RH","RH")
metric_labels = c("SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "SA","SA","SA","SA",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT",
                  "CT","CT","CT","CT")

df2plot$region = new_region_labels
df2plot$hemisphere = hemi_labels
df2plot$metric = metric_labels


colors2use = get_ggColorHue(7)

# SA F-stat
d2plot = subset(df2plot, df2plot$metric=="SA" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))

# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot, 
         display_numbers = round(mat2plot,digits=2), 
         number_color = "black", fontsize_number = 18,
         show_rownames=TRUE,
         labels_col = colnames(mat2plot),
         color = colorRampPalette(c('light blue','white','red'))(100),
         cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_manual(values = colors2use[1:3]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_manual(values = colors2use[4:7]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
                c("#785e00","#c49a00","#ffcc12"),
                c("#306800","#53b400","#76ff01"),
                c("#007459","#00c094","#0effc8"),
                c("#007b9f","#00b6eb","#39d2ff"),
                c("#6a3eff","#a58aff","#e0d7ff"),
                c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(0.6,1.3), c(0.8,1.15), c(0.75,1.35), c(0.8,1.35),
             c(0.8,1.4), c(0.8, 1.15), c(0.7,1.1))
                
for (ireg in 1:length(sig_reg)){
  region = sig_reg[ireg]
  reg_name = reg_names[ireg]
  metric_name = metric_names[ireg]
  p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
  p1 = p1 + geom_jitter(size = dotSize) +
    geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
    scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) + 
    xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) + 
    theme(text = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
  print(p1)
  
  tmp_df = d2plot[region,]
  if (tmp_df$hemisphere=="LH"){
    HEMI = "left"
  } else if (tmp_df$hemisphere=="RH"){
    HEMI = "right"
  }
  p = ggseg(.data = tmp_df, 
            atlas = chenAr, 
            mapping=aes(fill=label2fill), 
            position="stacked", 
            colour="black",
            hemisphere=HEMI) +
    scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) + 
    theme(text = element_text(size=fontSize))
  print(p)
}

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "CortexVol"
reg_name = "Total Cortical Volume"
metric_name = "Volume"
ylims = c(450000,700000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
p1

cols2use = c("#3d3d3d","#828282","#b3b3b3")
region = "SAtotal"
reg_name = "Total Surface Area"
metric_name = "SA"
ylims = c(90000,220000)
p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, fill = "Group"))
p1 = p1 + geom_jitter(size = dotSize, colour="black", pch=21) +
  geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
  scale_fill_manual(values=cols2use) + guides(colour = FALSE, fill = FALSE) + 
  xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims) + 
  theme(text = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize),
      plot.title = element_text(size=fontSize,hjust=0.5))
p1

# SA es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# SA es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# SA es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="SA" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenAr, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT F-stat
d2plot = subset(df2plot, df2plot$metric=="CT" & df2plot$fdrq<=fdr_thresh)
d2plot$label2fill = factor(c(1:dim(d2plot)[1]))

# make figure
mat2plot = t(d2plot[,c("es.TD_vs_Good","es.TD_vs_Poor","es.Good_vs_Poor")])
pheatmap(mat2plot, 
         display_numbers = round(mat2plot,digits=2), 
         number_color = "black", fontsize_number = 18,
         show_rownames=TRUE,
         labels_col = colnames(mat2plot),
         color = colorRampPalette(c('light blue','white','red'))(100),
         cluster_rows = TRUE, cluster_cols = TRUE)

# LH
df2use = subset(d2plot, d2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_manual(values = colors2use[1:4]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(d2plot, d2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          # mapping=aes(fill=Fstat), 
          mapping=aes(fill=label2fill), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_manual(values = colors2use[5:7]) + guides(fill=FALSE) + 
  # scale_fill_gradientn(colors = colorRampPalette(c("white","red"))(100), limits = c(0,15)) +
  # guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

sig_reg = rownames(d2plot)
reg_names = d2plot$region
metric_names = d2plot$metric
cols2use = list(c("#f53224","#f8766d","#fcbbb6"),
                c("#785e00","#c49a00","#ffcc12"),
                c("#306800","#53b400","#76ff01"),
                c("#007459","#00c094","#0effc8"),
                c("#007b9f","#00b6eb","#39d2ff"),
                c("#6a3eff","#a58aff","#e0d7ff"),
                c("#f916c4","#fb61d7","#fdacea"))
ylims = list(c(-0.4,0.1), c(-0.4,0.4), c(-0.4,0.3), c(0,0.7), 
             c(-0.3,0.3), c(-0.2,0.5), c(0,0.8))
                
for (ireg in 1:length(sig_reg)){
  region = sig_reg[ireg]
  reg_name = reg_names[ireg]
  metric_name = metric_names[ireg]
  p1 = ggplot(data = fs_data, aes_string(x = "Group", y = region, colour = "Group"))
  p1 = p1 + geom_jitter(size = dotSize) +
    geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) +
    scale_colour_manual(values=cols2use[[ireg]]) + guides(colour = FALSE) + 
    xlab("Group") + ylab(metric_name) + ggtitle(reg_name) + ylim(ylims[[ireg]]) + 
    theme(text = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize),
        plot.title = element_text(size=fontSize,hjust=0.5))
  print(p1)
  
  tmp_df = d2plot[region,]
  if (tmp_df$hemisphere=="LH"){
    HEMI = "left"
  } else if (tmp_df$hemisphere=="RH"){
    HEMI = "right"
  }
  p = ggseg(.data = tmp_df, 
            atlas = chenTh, 
            mapping=aes(fill=label2fill), 
            position="stacked", 
            colour="black",
            hemisphere=HEMI) +
    scale_fill_manual(values = colors2use[ireg]) + guides(fill=FALSE) + 
    theme(text = element_text(size=fontSize))
  print(p)
}

# CT es.Good_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.Good_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT es.TD_vs_Good
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Good), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf

# CT es.TD_vs_Poor
# LH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="LH")
pl = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="left") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) + 
  theme(text = element_text(size=fontSize))

# RH
df2use = subset(df2plot, df2plot$metric=="CT" & df2plot$hemisphere=="RH")

pr = ggseg(.data = df2use, 
          atlas = chenTh, 
          mapping=aes(fill=es.TD_vs_Poor), 
          position="stacked", 
          colour="black",
          hemisphere="right") +
  scale_fill_gradientn(colors = colorRampPalette(c("blue","white","red"))(100), limits=c(-0.8,0.8)) +
  guides(fill = guide_colourbar(nbin = 100)) +
  theme(text = element_text(size=fontSize))

pf = (pr)/(pl)
pf